基于R语言的微生物群落组成多样性分析——RDA分析

此前文章中小编已经给大家介绍了微生物群落组成多样性分析的三大主要板块——Alpha多样性分析、Beta多样性分析及群落结构组成分析。然而，在微生物高通量测序数据分析过程中，分析环境因子是否对微生物的群落结构造成一定影响也是其中很重要的一个板块。所以，小编今天要给大家讲的内容就是使用RDA分析方法分析环境因子对微生物群落的影响。 RDA，即冗余分析（Redundancy Analysis），是一种响应变量矩阵与解释变量之间多元多重线性回归的拟合值矩阵的PCA分析，说的通俗点其实就是约束化的PCA分析。通过RDA分析可以实现将样本和与其对应的环境因素同时反映在同一个图中，以此揭示出环境因子是否对样本分布或则样本微生物群落结构等造成影响。

加载包

rm(list=ls())#clear Global Environmentsetwd('D:\\桌面\\RDA')#设置工作路径
#加载包library(vegan)library(ggplot2)

加载数据

#OTU表格df <- read.table("otu.txt",sep="\t",header = T,row.names = 1,check.names = F)head(df)#环境因子数据env <- read.table("env.txt",sep="\t",header = T,row.names = 1,check.names = F)head(env)

RDA分析

1）使用decorana函数检测我们的数据是否符合做RDA的要求

print(decorana(t(df)))#根据DCA1的Axis Lengths值进行选择，如果>4.0选CCA；如果在3.0-4.0之间，选RDA和CCA都可以；如果<3.0, 选择RDA分析即可。

2）vegan包rda函数支持RDA分析

RDA <- rda(t(df),env,scale = T)

3）提取数据

#提取数据df_rda <- data.frame(RDA$CCA$u[,1:2],rownames(env))colnames(df_rda)=c("RDA1","RDA2","samples")# 提取物种得分df_rda_score <- data.frame(RDA$CCA$v[,1:2])#计算轴标签数据（=轴特征值/sum(所有轴的特征值)）RDA1 =round(RDA$CCA$eig[1]/sum(RDA$CCA$eig)*100,2)RDA2 =round(RDA$CCA$eig[2]/sum(RDA$CCA$eig)*100,2)

绘制RDA图

#读入分组文件group <- read.table("group.txt", sep='\t', header=T)#修改列名colnames(group) <- c("samples","group")#将绘图数据和分组合并df_rda <- merge(df_rda,group,by="samples")color=c("#1597A5","#FFC24B","#FEB3AE")#颜色变量p1<-ggplot(data=df_rda,aes(x=RDA1,y=RDA2, color=group))+#指定数据、X轴、Y轴，颜色 theme_bw()+#主题设置 geom_point(size=3,shape=16)+#绘制点图并设定大小 theme(panel.grid = element_blank())+ geom_vline(xintercept = 0,lty="dashed")+ geom_hline(yintercept = 0,lty="dashed")+#图中虚线 geom_text(aes(label=samples, y=RDA2+0.03,x=RDA1+0.03, vjust=0),size=3)+#添加数据点的标签 # guides(color=guide_legend(title=NULL))+#去除图例标题 labs(x=paste0("RDA1 (",RDA1,"%)"), y=paste0("RDA2 (",RDA2,"%)"))+#将x、y轴标题改为贡献度 stat_ellipse(data=df_rda, level=0.95, linetype = 2,size=0.8, show.legend = T)+ scale_color_manual(values = color) +#点的颜色设置 scale_fill_manual(values = c("#1597A5","#FFC24B","#FEB3AE"))+ theme(axis.title.x=element_text(size=12),#修改X轴标题文本 axis.title.y=element_text(size=12,angle=90),#修改y轴标题文本 axis.text.y=element_text(size=10),#修改x轴刻度标签文本 axis.text.x=element_text(size=10),#修改y轴刻度标签文本 panel.grid=element_blank())#隐藏网格线p1

添加环境因子数据

#提取环境因子得分df_rda_env <- RDA$CCA$biplot[,1:2]df_rda_env <- as.data.frame(df_rda_env)head(df_rda_env)# 添加环境因子数据p1+ geom_segment(data=df_rda_env,aes(x=0,y=0,xend=df_rda_env[,1],yend=df_rda_env[,2]), color="black",size=0.8, arrow=arrow(angle = 35,length=unit(0.3,"cm")))+ geom_text(data=df_rda_env,aes(x=df_rda_env[,1],y=df_rda_env[,2], label=rownames(df_rda_env)),size=3.5, color="blue", hjust="inward", vjust=0.5*(1-sign(df_rda_env[,2])))+  theme(legend.position = "top")

AI美化

源码及作图数据在后台回复“RDA”即可获取！！！

爱我请给我好看!